3. 1群の正規分布の分析

https://gyazo.com/4f7dfd91756d579cff6ea820ee6ff329

https://amzn.to/32dVAdF

3.1 事後分布の近似

前章で学んだように、統計的推測の対象となる母数の事後分布は、一般的に式表現することは容易

しかし、平均や標準偏差を求めるなど、性質を評価できる形式で事後分布を表現することはほとんどの場合に不可能

事後分布を評価できないという理由から、ベイズ統計学が20世紀中に市民権を得ることはできなかった

近年、計算機パワーを利用し、発想の転換を図ることで、事後分布の詳細な性質を評価する事が可能になった

3.1.1. MCMC法

事後分布に従う母数を乱数として発生させ、事後分布をあたかもデータ分布のように入手すること

マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo method, MCMC法)

本書ではMCMC法の一つの手法であるハミルトニアンモンテカルロ法(Hamiltonian Monte Carlo method, HMC法)を利用する

HMC法は、物理学分野の力学的エネルギーの原理を応用した方法

現代的なベイズ分析では、統計的推論の多くをMCMC法でサンプリングされた母数の分析を通じて行う

3.1.2. 生成・チェイン・バーンイン

MCMC方は、同時事後分布に従う乱数を、蛇口から水が流れるように経時的に生成する

第$ m期に発生した乱数を$ \bm \theta^{(m)}(m = 1, \cdots , M)と表記する

バーンイン(burn-in, 焼き入れ)期間、ウォームアップ(warmup)期間

ただし第$ 1期から第$ M期まで生成させても、初期の乱数は同時事後分布に従わない

そこで最初から$ B期までの乱数を捨てて利用しないこととする

事後分布の性質を調べるためには、バーンイン以後($ m=B+1, B+2, \cdots, M)の有効な乱数を用いる

バーンイン期間($ m = 1, 2, \cdots, B)の乱数は捨て、以後一切の分析に使用しない

チェイン(chain)

乱数列

乱数列の数をチェイン数という

$ M = 21000のチェインを5つ発生させ、バーンイン期間を$ B = 1000とし、得られた$ T= 100,000($ =(21000-1000) \times 5)個の有効な乱数を用いて描いた散布図が図2-1だった

これを1列に並べて$ \bm \theta^{(t)}(t = 1, \cdots, T)と表記する

3.1.3. トレースプロット

トレースプロット(trace plot)

経時$ mに沿って乱数の値を折れ線で表現したグラフ

図2-1の事後分布を求めたときのトレースプロット

https://gyazo.com/6a5991e495f2ee78f1ce959717d2bb39

事後分布から乱数が発生している必要条件

もはや折れ線には見えない、真横に広げたテープのような形状が観察されること

「トレースプロットがドリフト(drift)している」

https://gyazo.com/d643963185cdf0b36ab3a6c42001a166

登ったり降りたりの形状

目的とする確率分布に従った乱数を得ていない視覚的特徴

同時に描かれているチェインの折れ線が互いに区別がつかないことも、事後分布から乱数が発生している必要条件

上部と下部のように、いくつかの折れ線が違った水平レベルで描かれているなら、それぞれの乱数が目的とする確率分布をカバーしていない視覚的証拠

https://gyazo.com/372cd508a568ba2ab7a82ca2ef7c8bab

このようにトレースプロットは「目的としている事後分布から乱数が正しく発生しているか否か」の視覚的評価に利用できる

3.1.4. 収束判定指標

トレースプロットの観察から正常と考えられる倍位には、乱数列の数値的評価の指標を参照する

https://gyazo.com/c064ccb27f6b6806472eb1ec41119efc

ここで示されているのは収束判定指標$ \hat{R}と有効標本数$ n_{eff}

$ x^*に関しては後述

収束判定指標$ \hat{R}

チェインの間の散らばりがチェイン内の散らばりに比べて大きい場合には、事後分布から正しく乱数が発生していないことが疑われ、そのことを$ \hat Rが警告してくれる

チェイン数が$ 1の場合には、それを複数のチェインに分割することで計算する

$ \hat Rは$ 1.1ないし$ 1.2以下であればよいとされる

3.1.5. 有効標本数

有効標本数(effective sample size)$ n_{eff}

確率分布を近似する乱数は、互いに関係し合わないことが理想的

しかし、MCMC法は経時的に乱数を発生するから、現実的には$ m番目と$ m+1番目の乱数は関係を持っている

有効標本数は「生成された乱数が理想的に無関係である乱数の何個分に相当するか」の推定値

ここでは事後分布を近似するために、$ 100000個の乱数を残した。

その$ \muに関する乱数は、理想的な独立した乱数の$ 55404個分に相当すると推定されている

その$ \sigmaに関する乱数は、理想的に独立した乱数の$ 55882個分に相当すると推定されている

3.1.6. 推定値の不安定性

MCMC法の欠点は、事後分布の評価が完全には定まらないこと

MCMCの手法を変更したり、アルゴリズムを変えたり、ソフトウェアを更新すればもちろんのこと、乱数の種を変えるだけでも、事後分布の%点が変化する

平均や標準偏差の推定値も変化する

本書に登場する統計モデルは単純なので、分析結果の解釈が影響を受けるほどの不安定さにはならないことがほとんど

しかしより安定させたい場合には、大きな$ Mを指定し、結果として$ Tを大きくする

3.2 推定量・精度

データに関する情報はデータ分布に全て含まれていた

同様に母数に関する情報は事後分布にすべて含まれている

データ分布と事後分布を対応づけながら、母数の推測に入門する

3.2.1. 点推定量

データ分布を勉強したときには、要約統計量を求め、データの有する情報を縮約的に記述した

母数に関しても、常に事後分布を参照するのは煩雑なので同じように縮約する

データ分布では、まず位置を表現する代表値を計算した

平均値$ \bar x・中央値$ x_{med}・最頻値$ x_{mod}

つまり、データ分布を点で代表させた

同様に母数の事後分布に関しても、まずは点で代表させる

点推定(point estimation)

母数を点で代表させること

点推定量(point estimator)

推定の方法

点推定値(point estimate)

推定された値

母数$ \thetaの事後分布も、データの分布と同様に、平均値・中央値・最頻値という3つの観点から代表させる

それらが以下に示す代表的な母数の点推定量になる

事後期待値(expected a posteriori, EAP $ \theta_{eap})

事後分布の平均値であり、$ \theta^{(t)}の平均値を推定値とする

$ \hat \theta_{eap} = \frac{1}{T}(\theta^{(1)} + \theta^{(2)} + \cdots + \theta^{(T-1)} + \theta^{(T)})

この方法をEAP推定量という。

ただし、平均値の計算には、バーンイン以後の乱数しか使用していないことに留意

平均値ばかりでなく以後のすべての分析・作図においてもバーンイン以後の乱数しか使用しない

事後中央値(posterior median, MED, $ \theta_{med})

事後分布の中央値

$ \theta^{(t)}を小さい順にソートして、その中央の値をMED推定量とする

$ \hat \theta_{med} = \begin{cases} \frac{T+1}{2}番目の乱数 & Tが奇数の場合 \\ \frac{T}{2}番目と\frac{T}{2} + 1番目の乱数の平均 & Tが偶数の場合 \end{cases}

要するに50%点である

事後確率最大値(maximum a posteriori, MAP, $ \theta_{map})

事後分布のモード

$ \theta^{(t)}のヒストグラムの最も度数の大きい階級の階級値を推定値に使うことをMAP推定量($ \hat\theta_{map})という

3.2.2. 事後標準偏差

データ分布を学習した際には、代表値を観察した後、分散$ s^2・標準偏差$ sという2つの要約統計量で散布度を評価した

正規分布の散布度は母分散$ \sigma^2と母標準偏差$ \sigmaだった

事後分布の散布度の小ささは点推定値の精度である

事後分布の散布度は事後分散と事後標準偏差で評価する

事後分散(posterior variance, $ \hat\sigma^2_\theta)

事後分布の分散

$ \hat \sigma^2_\theta = \frac{1}{T}((\theta^{(1)}-\hat\theta_{eap})^2 + \cdots (\theta^{(t)}-\hat\theta_{eap})^2) + \cdots + (\theta^{(T)}-\hat\theta_{eap})^2

事後標準偏差(posterior standard deviation, post.sd, $ \sigma_\theta)

事後分布の標準偏差

事後分散の推定値の平方根で推定する

$ \hat\sigma_\theta = \sqrt{\hat\sigma^2_\theta}

事後標準偏差は$ \thetaの標準偏差である

その推定値$ \hat\sigma_\thetaは、母数$ \thetaが$ \hat\theta_{eap}の周辺でどれほど散らばっているかの精度の指標として利用する

EAP推定値のまわりに事後分布が密集していれば、母数の評価値としてEAP推定値を利用しても不都合は生じにくいだろう

しかし事後分布が広範囲に渡っている場合には、EAP推定値の代表性が疑われる

母数の推定が安定的に行われているとはいえない

事後標準偏差が大きいなら、$ nが小さいということ

事後標準偏差を小さくするには、データを追加する必要がある

3.2.3. 確信区間

区間推定(interval estimation)

特定の値で母数を評価する点推定に対して、幅をもたせて区間で母数の評価を行う方法

事後分布の両端から$ \frac{\alpha}{2}\%の面積を切り取って残った中央部の$ (1-\alpha)\%の面積に対応する区間を$ (1-\alpha)\%の両側確信区間(credible interval)という

下側(上側)から$ \alpha\%の面積を切り取って残った$ (1-\alpha)\%の面積に対応する区間を$ (1-\alpha)\%上側(下側)確信区間という

上側・下側確信区間をまとめて片側確信区間という

本書では$ \alpha = 5\%として$ 95\%確信区間を利用する

確信区間は信用区間と呼ばれることもある

ベイズ統計学では確信区間ばかりでなく、最高事後密度区間(highest posterior density interval, HPDI)を利用することもある

$ A\%最高事後密度区間は、事後分布の密度が高い部分の$ A\%の範囲と定義される

確率の大きい順に、合計$ A\%になるまで階級を区間に組み込んだもの

本書では割愛

3.3 予測分布

手元のデータではなく、将来観測されるであろうデータ$ x^*を予測したい場合がある

具体的に「知覚時間」で例をあげるならば、実験が終わった後に、同じ手続きで同じ人がもう一度トライした場合の知覚時間

予測分布(predictive distribution)

将来観測されるであろうデータ$ x^*の分布

予測分布には2種類ある

3.3.1 事後予測分布

事後予測分布(posterior predictive distribution)$ f(x^*|\bm x)

観測してしまったデータ$ xを所与とした時の未来のデータ$ x^*の条件付き分布であり、以下の乱数列で近似する

$ x^{*(t)} \sim f(\theta^{(t)})

正規分布モデル

$ x^{*(t)} \sim N(\mu^{(t)}, \sigma^{(t)})

データ生成分布は正規分布であったが、母数$ \muと$ \sigma自体が事後に分布し、確率的に揺れるから、事後予測分布は必ずしも正規分布にならない

事後分布の頻度(確率密度)に応じて母数$ \muと$ \sigmaを呼び出し、それを利用して乱数を発生させる

これが将来のデータを予測するための事後予測分布

長所

事後分布の情報をあまさず利用する精密さ

短所

単純なデータ生成分布として表現できないこと

事前データ・事後分布を持ち続ける必要があり、これは煩雑

事後予測分布は、将来を丁寧に予測したい場合に利用するとよい

3.3.2. 条件付き予測分布

条件付き予測分布(conditional predictive distribution)

モデルの分布を$ f(x|\bm \theta)とした場合の条件付き予測分布

$ f(x^*|\hat\bm\theta)

ここで$ \hat\bm\thetaには何らかの点推定値を用いる

条件付き予測分布は、母数の推定値$ \hat\bm\thetaを所与とした時の未来のデータ$ x^*の条件付き分布である

条件付予測分布の(メタ分布)を求める時には点推定値の代わりに$ \bm\theta^{(t)}を用いる

長所

将来の予測が点推定値にだけ依存することによる取り扱いの容易さ

条件付き予測分布はデータ生成分布(この場合は正規分布)になるから扱いやすい

短所

事後分布の豊かな情報を点推定値だけで要約していることによる情報損失

3.3.3. 予測区間

予測分布の両端から$ \frac{\alpha}{2}\%の面積を切り取って残った中央部の$ (1-\alpha)\%の面積に対応する区間を$ (1-\alpha)\%両側予測区間(prediction interval)という

予測分布の下側(上側)から$ \alpha\%の面積を切り取って残った$ (1-\alpha)\%の面積に対応する区間を$ (1-\alpha)\%上側(下側)予測区間という

上側・下側予測区間をまとめて片側予測区間という

3.4. 「知覚時間」のベイズ的推測

「知覚時間」の母数に関するベイズ的推測を行う

ベイズ統計学に限らず、統計的推測をする際には研究目的(research object)、あるいは研究上の問い(research question, リサーチ・クエスチョン), RQを重視し、それを自覚することが大切

例えば以下のようなRQが着想される

RQ.1 平均値の点推定

e.g. 「知覚時間」の平均$ \muは何秒だろう

RQ.2 平均値の両側区間推定

e.g. 平均$ \muはどの区間にあるだろう

RQ.3 平均値の片側区間推定

e.g. 平均$ \muは高々あるいは少なくとも何秒だろう

RQ.4 標準偏差の点推定・区間推定

e.g. 平均的な知覚の散らばり$ \sigmaは何秒だろう。それはどの区間にあるだろう

RQ.5 予測区間

e.g. 同じ手続きで同じ人がもう一度実験したら、$ 95\%の確率で、何秒から何秒の間に測定値は収まるだろう

3.4.1. 平均の推測

母数$ \muの事後分布

https://gyazo.com/48d993d3c76f2611a7d657a0961d9837

中央の白い棒と灰色の棒を合わせた区間が$ 95\%両側確信区間

下の灰色の棒と黒い棒を除いた区間が$ 95\%上側確信区間

上の灰色の棒と黒い棒を除いた区間が$ 95\%下側確信区間

母数の事後分布と事後予測分布の要約統計量

https://gyazo.com/011de60e97160082d1f52560203beeeb

$ \muに関するEAP、MEDは$ \hat \mu_{eap} = 31.04秒であり、$ \hat \mu_{med} = 31.04用である(RQ.1への回答)

事後分布の積率系の統計量と分位系の統計量からだけではMAP推定値は求まらない

ただし、$ f(\bm x|\mu, \sigma)(1/(\beta_\mu - \alpha_\mu))(1/(\beta_\sigma - \alpha_\sigma)) \propto f(\bm x|\mu, \sigma)によりMAP推定値がMLEに一致する

さらにMLEはデータの標本平均に一致することが知られているから$ \bar x = \hat \mu_{mle} = \hat \mu_{map} = 31.04秒である

少数第2位までは、EAP、MED、MAPに差はなかった

母平均の点推定値は$ 31.04秒であり、その平均的な推定誤差であるpost.sdは$ \hat \sigma_\mu = 0.52秒であった

post.sdの小ささから、$ \muは$ 30秒より上にあると推定される

$ \muの確信区間を求めてこの問いを精査する

$ \muを初めとする一般の母数の事後分布は正規分布するとは限らないから、確信区間を求めるためには、それに対応する$ \theta^{(t)}の%点を参照する火チウ用がある

表3-2より、2.5%点と97.5%点を参照し

$ p(2.5\%点 \leq \mu \leq 97.5\%点) = 0.95

であるから、$ \mu の$ 95\% 両側確信区間は$ [30.02秒, 32.06秒] であることがわかる

母平均$ \muは$ 95\%の確率で固定されたこの区間に入っていると解釈する(RQ.2への回答)

100分の2秒であるが、30秒からは外れている

表3-2より「研究仮説(research hypothesis)あるいは命題(proposition): $ \mu < 31.88秒」が真である確率は95%である(RQ.3 への回答)

同様に「研究仮説: $ 30.20秒 < \mu」が真である確率は$ 95\%である

表3-2の下から$ 5\%点は、不等号の向きを変えれば上から$ 95\%として解釈できる

これが$ \muの$ 95\%片側確信区間から導かれる知見

3.4.2. 標準偏差の推測

$ \sigmaに関するEAP, MEDは表3-2より$ \hat\sigma_{eap}=2.27秒、$ \hat\sigma_{med}=2.22秒である(RQ.4の前半への回答)

MAPは表からはわからないが、$ f(\bm x|\mu, \sigma)(1/(\beta_\mu - \alpha_\mu))(1/(\beta_\sigma - \alpha_\sigma)) \propto f(\bm x|\mu, \sigma)によりMLEに一致する

さらにMLEはデータの標準偏差に一致することが知られているから$ s=\hat\sigma_{mle}=\hat\sigma_{map}=2.07秒

図3-4は$ \sigmaの事後分布

https://gyazo.com/12b63949406b06e6e9f376ceb8b7b53b

平均$ \muとは異なり、標準偏差$ \sigmaは推定量によって違いがある

標本標準偏差とEAPは一般的に値が異なる

これは$ \sigmaの事後分布が正に歪んでいるため

まれに標準偏差はとても大きな値をとる

このため標本標準偏差に一致する分布の頂点ではなく、それより大きい側に寄った点を敢えて推定値とするEAPの方針はその意味で合理的

事後分布に限らず、分布が非対称で右に長く裾を引いている形状のとき、その分布は正に歪んでいると表現する

事後分布が正に歪んでいる場合は、一般的にMAP<MED<EAPとなる

図3-4は右に裾を引いているが、これは「知覚時間」に限ったことではない

標準偏差の事後分布は、常に正に歪んだ形状となり、MAP<MED<EAP

また「知覚時間」における$ \sigma の$ 95\% 両側確信区間は$ [1.65秒, 3.22秒] である(RQ.4の後半への回答)

3.4.3. 事後予測分布

https://gyazo.com/2853242abb62fdd43c3aae03c37f1ccf

中央の白い棒と灰色の棒を合わせた区間が$ 95\%両側予測区間

下の灰色の棒と黒い棒を除いた区間が$ 95\%上側予測区間

上の灰色の棒と黒い棒を除いた区間が$ 95\%下側予測区間

実験が終わった後に、同じ手続きで同じ人がもう一度トライした知覚時間の目安として、$ x^*の事後予測分布の$ 95\%両側予測区間を利用できる

それは表3-2より、$ [26.35秒, 35.71秒]である(RQ.5への回答)

$ \mu の$ 95\% 確信区間$ [30.02秒, 35.71秒] よりもだいぶ広くなった

$ x^*の$ 95\%予測区間が$ \muの確信区間よりも広くなることは一般的性質である

この場合、$ x^*は将来の測定値であるから、母平均$ \muよりもさらに長くなることも、さらに短くなることもある

むしろ注意すべきことは、標本平均と標本標準偏差を用いて1. データ分布の要約で導いた$ F(\mu + 1.96\sigma|\mu, \sigma) - F(\mu-1.96\sigma|\mu, \sigma) \simeq 0.95 式の$ 95\% 予測区間$ [25.99秒, 35.10秒] よりも、事後予測区間の方が広いということ

これは$ \muや$ \sigmaが分布しているためである

言い換えるならば、母数の推定に伴う不確実性が加味されるために、事後予測分布による予測区間の方が広くなる

3.4.4. 条件付き予測分布

条件付き予測分布は、EAP, MED, MAPの各推定値を利用し、正規分布の密度関数を用いてそれぞれ以下のようになる

$ f(x^*|\hat\mu_{eap} = 31.04, \hat\sigma_{eap} =2.27)

$ f(x^*|\hat\mu_{med} = 31.04, \hat\sigma_{med} = 2.22)

$ f(x^*|\hat\mu_{map} = \hat\mu_{mle} = 31.04, \hat\sigma_{map} = \hat\sigma_{mle} = 2.07)

標本平均と標本標準偏差を利用した正規分布へのあてはめは、MAP推定値による条件付き予測分布に一致する

「知覚時間」の条件付き予測分布は、正規分布であるから図示は割愛する

→4. 生成量と研究仮説が正しい確率